#ifndef _QGCM_PSYSTEM_
#define _QGCM_PSYSTEM_
#include "datatypes.h"
#include "constants.h"
#include "forward_classes.h"
#include <mpParser.h>
#include <netcdfcpp.h>

using namespace mup;

class psystem  {
	friend class model_class;
	friend class ssystem;
	friend class compiler;
	friend class debugger;

	private:
	struct {
		bool time_initialized;
		bool IC_initialized;
	} model_state;


	int numfields;
	psystem(const psystem &); // no copy constructor should be callled!
	psystem & operator=(const psystem &); // no equal assignments allowed

	model_class * mo; // pointer to owner class

	vector<string> variable_names; 
	vector<Variable> variables;
	vector< vector<Value> > stored_variables; // keep old values in case we need to rewind

	vector<string> function_names;
	vector<ParserX*> functions; // functions are evaluated in the order they appear here
	vector<Variable> function_values; // values are stores in this array at each time step
	// . . . . . . . . . . .
	vector<string> sfield_names; // names of fields evaluated by the s-subsystem
	vector<Variable> sfields;  // evaluated by s-system.
	vector< rank2 > sfield_indices; // field indices corresponding to the flux element
	// * * * * * * * * 
	vector<ParserX * > mean_field_initializers; // used to create initial conditions
	vector<Variable> mean_fields; // stores mean_fields themselves
	vector<Value> mean_field_temp;
	vector< vector<Value> > stored_meanfields;
	vector<Variable> mean_field_derivatives; // stores pointers to output of functions that
	// are derivatives of the respective fields 
	vector<Value> mean_field_deriv_temp; // temporary storage
	// vector to be evaluated each time
	vector<string> mean_field_names;

	int tindex;
	unsigned int numtsteps; 


	private:
	void save_state(int tstep);

	void Eval(int ti, bool provisional);

	bool restore_model_state(NcFile * ncfile);

	public:
	psystem();
	~psystem();
	void init(model_class * mo);
	void init_state(); // set initial conditions
	void addfluxfield(const string &, INDEX, INDEX, const int );

	const int get_tindex() const;
	void init_tindex(const unsigned int &);

	void rewindtime(const unsigned int &);
	void adjusttime(const unsigned int &); // adjust to specific time-step

	bool save_model_state(NcFile *) const;

	const unsigned int getvariablesize() const;
	const Variable & getvariablebyindex(const unsigned int) const;
	const string & getvariablename(const unsigned int) const;
	bool doesvariableexist(const string & name) const;
	const vector< vector< Value > > & getvariablehistory() const;


	const Variable & getmeanfield(unsigned int index) const;
	const unsigned int getmeanfieldsize() const;
	const Variable & getmeanfieldderivative(unsigned int index) const;
	const string & getmeanfieldname(unsigned int index) const;
	const string & getmeanfieldinitexpr(unsigned int index) const;
	const vector< vector< Value > > & getmeanfieldhistory() const;

	int getfieldbyname(const string &) const;

	void assignparserdata(ParserX &) const;

	const Variable & getsmeanfield(unsigned int index) const;
	const unsigned int getsmeanfieldsize() const;
	const string & getsmeanfieldname(unsigned int index) const;
	int getsfieldbyname(const string &) const;
	const rank2 & getsmeanfieldindices(unsigned int index) const;

	const Variable & getfunctionvalue(unsigned int index) const;
	const unsigned int getfunctionssize() const;
	const string & getfunctionname(unsigned int index) const;
	const string & getfunctionexpr(unsigned int index) const;
	int getfuncbyname(const string &) const;


	void stepforward(REAL timestep, bool nosys, bool nonlinear);
	void psystem_stepforward_init(Variable * & ,Variable *  &);

	int getnumtsteps() const;
};


#endif
